function mat = Spin( j )
% calculation of the spin matrix with spin angular momentum j (3D double array)
% p : raising operator \J_{+}
% using the formula \J_{+} |j,m> = sqrt(j(j+1)-m(m+1))
%SPIN Summary of this function goes here
%   Detailed explanation goes here

    p=diag(sqrt((1:2*j).*(2*j:-1:1)),1);
    mat(:,:,1)=(p+p')/2;
    mat(:,:,2)=(p-p')/(2*1i);
    mat(:,:,3)=diag(j:-1:-j);
end

